\documentclass[clear=false]{groupreport}


\special{dvipdfmx:config z 0}
% % 取消pdf的压缩，最终版取消这个命令即可

% \title{柱壳的振动特性分析}%标题
% \author{BY2104502 杜晨鸿}%作者
% \date{2021-12}%日期
 
%——————以下宏包功能不一定用上，请自行学习——————————
% \usepackage{natbib}
% \usepackage[table,xcdraw]{xcolor}%提供更多的色彩
\usepackage{extarrows} %长等号扩展
\usepackage{booktabs}%三线表格
\usepackage{amsmath}%扩充数学富符号
\usepackage{multirow}
% \usepackage{cite}
\usepackage{graphicx}%引入图片宏包
% \usepackage{subfigure}%小图片
\usepackage{listings}  %代码块当然数学作业也用不到
\usepackage{makecell}


%不熟悉latex的数学公式代码可以去www.latexlive.com 生成。
% 设置题注格式


\usepackage{chngcntr}
\counterwithout{equation}{section}
\counterwithout{figure}{section}

\usepackage{bm}

% \renewcommand{\cleardoublepage}{}
% \renewcommand{\clearpage}{}

\begin{document}%正文部分

% 无法加在cls里的全局设置
\setmainfont{times.ttf}[
  ItalicFont = timesi.ttf,
  BoldFont = timesbd.ttf,
  BoldItalicFont =timesbi.ttf
]
\zihao{-4}

\frontmatter

\captionsetup{labelformat=default,labelsep=space}  % 去除题注冒号

\includepdf{cover.pdf}
% \newpage%新的一页
\thispagestyle{empty}
\leavevmode\newpage
\thispagestyle{empty}
% \includepdf{sign.pdf}
% \thispagestyle{empty}
% \leavevmode\newpage
% \thispagestyle{empty}

\pagenumbering{gobble}
\tableofcontents%生成目录
\clearpage
% \pagestyle{nofoot}
\newpage%新的一页

\mainmatter{
\chapter{选定几何及材料参数}
% \thispagestyle{fancy}

\punctstyle{kaiming}
柱壳\songkuohao{示例括号}的几何及材料参数见表\ref{materal}。



\begin{table}[H]
  \centering
  \caption{几何及材料参数}\label{materal}
  % \begin{center}
  \begin{tabular}{lll}
    \toprule
    参数   & 符号     & 值                          \\\midrule
    壳厚   & $h$    & 1mm                        \\
    长度   & $L$    & 100mm                      \\
    半径   & $R$    & 50mm                       \\
    密度   & $\rho$ & $\rm 8\times10^{-9}t/mm^3$ \\
    杨氏模量 & $E$    & 200000MPa                  \\
    剪切模量 & $G$    & 76923MPa                   \\
    泊松比  & $\nu$  & 0.3                        \\
    \bottomrule
  \end{tabular}
\end{table}

正体希腊字母为$\alpha\upalpha\pi\uppi\varepsilon\upvarepsilon$

\begin{align}
  \sum_{\uppi}^{\upbeta}
\end{align}

斜体times为：

\emph{abcdefghijklmnopqrstuvwxyz}

{\setmainfont{timesi.ttf} abcdefghijklmnopqrstuvwxyz}

\begin{align}
  abcdefghijklmnopqrstuvwxyz \\
  xyz\boldsymbol{x}G\boldsymbol{G}\varPhi\bm{\varPhi}
\end{align}


\chapter{有限单元法求解}

\section{有限元模型}
\subsection{3级示例标题}
{\zihao{-4}本文}所采用\songkuohao{示例123abc括号}的有限元模型如图\ref{modal-fem}所示。采用ANSYS中的8节点四边形厚壳单元SHELL281。
共有单元31400个，节点94829个。对左侧一圈节点施加固定约束，采用分块兰索法进行模态分析。
\begin{figure}[H]
  \centering
  \includegraphics[height=150pt]{./images/modal-fem.png}
  \caption{123有限元模型}\label{modal-fem}
\end{figure}


\section{计算结果}
在有限单元法算得的模态中，每个共振频率都对应2个频率相同，振型相差一定角度\songkuohao{周向}的振型。
本文仅记录前三节互异模态\songkuohao{即频率和振型均不相同}，如图\ref{level1-fem}-图\ref{level3-fem}所示，前三阶振型分别为6节径、8节径、4节径。其频率见表\ref{freq-fem}。
\begin{figure}[H]
  \centering
  \subfigure[\zihao{5}子图1号]{
    \includegraphics[height=150pt]{./images/cyc000.jpg}
  }
  \subfigure[\zihao{5}2号子图]{
    \includegraphics[height=150pt]{./images/cyc001.jpg}
  }
  \caption{有限元法一阶模态}\label{level1-fem}
\end{figure}

\begin{figure}[H]
  \centering
  \includegraphics[height=150pt]{./images/cyc003.jpg}
  \includegraphics[height=150pt]{./images/cyc002.jpg}
  \caption{2阶有限元法模态}\label{level2-fem}
\end{figure}

\begin{figure}[H]
  \centering
  \includegraphics[height=150pt]{./images/cyc004.jpg}
  \includegraphics[height=150pt]{./images/cyc005.jpg}
  \caption{有限元法三阶模态}\label{level3-fem}
\end{figure}

\begin{table}[H]
  \centering
  \caption{有限元解的前三阶互异模态频率}\label{freq-fem}
  \begin{tabular}{ll}
    \toprule
    阶数 & 频率\songkuohao{Hz} \\\midrule
    1  & 1491.7 \\
    2  & 1646.4 \\
    3  & 2324.5 \\
    \bottomrule
  \end{tabular}
\end{table}

\begin{table}
  \centering
  \caption{A segment of the ``Wide-angle eye'' dataset with direction strategy labels}\label{tbl:table4}
  \begin{tabular}{c m{7cm} c}
    \toprule
    \bfseries Index & \multicolumn{1}{c|}{\bfseries 5$\times$5 ``Wide-angle Eye'' Data} & \bfseries Direcion Strategy Label \\
    \midrule
    1               &
    \begin{equation}
      \left(
      \begin{array}{ccccc}
        116.81 & 21.14 & 109.43 & 72.00 & 1.90\nonumber \\
        115.83 & 18.83 & 109.67 & 71.30 & 1.86\nonumber \\
        115.51 & 18.27 & 109.91 & 69.68 & 1.87\nonumber \\
        115.83 & 21.76 & 109.46 & 70.83 & 1.90\nonumber \\
        115.83 & 21.30 & 110.87 & 72.22 & 1.90\nonumber \\
      \end{array}
      \right)
    \end{equation}
                    & \bfseries-1                                                                                           \\
    \midrule
    \bfseries 2     &
    \begin{equation}
      \left(
      \begin{array}{ccccc}
        115.83 & 18.83 & 109.67 & 71.30 & 1.86\nonumber \\
        115.51 & 18.27 & 109.91 & 69.68 & 1.87\nonumber \\
        115.83 & 21.76 & 109.46 & 70.83 & 1.90\nonumber \\
        115.83 & 21.30 & 110.87 & 72.22 & 1.90\nonumber \\
        115.83 & 21.14 & 111.57 & 72.00 & 1.89\nonumber \\
      \end{array}
      \right)
    \end{equation}
                    & \bfseries-1                                                                                           \\
    \midrule
    \bfseries 3     &
    \begin{equation}
      \left(
      \begin{array}{ccccc}
        130.42 & 42.75 & 119.93 & 84.26 & 2.07\nonumber \\
        128.47 & 40.27 & 120.63 & 84.50 & 2.07\nonumber \\
        122.64 & 38.89 & 120.41 & 82.18 & 2.07\nonumber \\
        122.00 & 38.89 & 120.43 & 81.71 & 2.07\nonumber \\
        121.67 & 38.89 & 120.67 & 81.25 & 2.06\nonumber \\
      \end{array}
      \right)
    \end{equation}
                    & -1                                                                                                    \\
    \midrule
    % \bfseries 4     &
    % \begin{equation}
    %   \left(
    %   \begin{array}{ccccc}
    %     92.50 & 48.61 & 113.57 & 67.36 & 1.95\nonumber \\
    %     92.50 & 45.52 & 113.81 & 67.89 & 1.94\nonumber \\
    %     92.50 & 44.44 & 114.06 & 67.59 & 1.92\nonumber \\
    %     92.82 & 43.20 & 114.07 & 67.82 & 1.95\nonumber \\
    %     95.42 & 45.37 & 114.09 & 67.82 & 1.91\nonumber \\
    %   \end{array}
    %   \right)
    % \end{equation}
    %                 & \bfseries 0                                                                                           \\
    % \midrule
    \bfseries 5     &
    \begin{equation}
      \left(
      \begin{array}{ccccc}
        92.50 & 45.52 & 113.81 & 67.89 & 1.94\nonumber \\
        92.50 & 44.44 & 114.06 & 67.59 & 1.92\nonumber \\
        92.82 & 43.20 & 114.07 & 67.82 & 1.95\nonumber \\
        95.42 & 45.37 & 114.09 & 67.82 & 1.91\nonumber \\
        95.42 & 45.99 & 113.87 & 68.98 & 1.95\nonumber \\
      \end{array}
      \right)
    \end{equation}
                    & \bfseries 0                                                                                           \\
    \midrule
    \bfseries 6     &
    \begin{equation}
      \left(
      \begin{array}{ccccc}
        137.67 & 42.57 & 109.90 & 65.11 & 1.98\nonumber \\
        137.67 & 44.09 & 108.84 & 65.36 & 2.00\nonumber \\
        149.97 & 37.95 & 109.09 & 65.36 & 2.00\nonumber \\
        135.50 & 48.77 & 108.84 & 65.61 & 1.97\nonumber \\
        122.12 & 48.78 & 109.09 & 65.61 & 1.96\nonumber \\
      \end{array}
      \right)
    \end{equation}
                    & \bfseries +1                                                                                          \\
    \bottomrule
  \end{tabular}
\end{table}

\section{假设模态法求解}

\subsection{动能及势能}
任意连续体的动能可以表示为其中所有质点动能的积分，即\upcite{lee2015free}
\begin{equation}
  T=\frac{\mathrm{1}}{\mathrm{2}}\rho\int_0^L\int_0^{\rm 2\pi}\int_{-h/{\rm 2}}^{h/{\rm 2}}
  (\dot{\xi}^2+\dot{\eta}^2+\dot{\omega}^2)R{\rm d}z{\rm d}\varphi{\rm d}x
\end{equation}
其中$x$、$\varphi$、$z$ \textit{z}分别表示轴向、周向、径向坐标，其中$\varphi$的量纲为弧度，其余为长度。
 $\xi$ 、 $\eta$ 、 $\omega$分别表示轴向、周向、径向位移，量纲均为长度。
对于薄壳，可假径向各质点具有相同的速度，即位移仅为$x$、$\varphi$的函数，与$z$无关，
于是可对其内层积分进行计算，将上式化简为
\begin{equation}
  T=\frac{hR\rho}{2}\int_0^L\int_0^{\rm 2\pi}
  (\dot{\xi}^2+\dot{\eta}^2+\dot{\omega}^2){\rm d}\varphi{\rm d}x
\end{equation}

任意线弹性连续体的变能可以表示为
\begin{equation}
  U=\frac{R}{2}\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  (\sigma_x\epsilon_x+\sigma_{\varphi}\epsilon_{\varphi}+\sigma_z\epsilon_z
  +\tau_{xz}\gamma_{xz}+\tau_{z\varphi}\gamma_{z\varphi}+\tau_{x\varphi}\gamma_{x\varphi})
  {\rm d}z{\rm d}\varphi{\rm d}x
\end{equation}
我们将其假设为平面应力问题，忽略径向的正应力及切应力，可得\upcite{lee2015free}
\begin{equation}
  U=\frac{R}{2}\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  (\sigma_x\epsilon_x+\sigma_{\varphi}\epsilon_{\varphi}
  +\tau_{x\varphi}\gamma_{x\varphi})
  {\rm d}z{\rm d}\varphi{\rm d}x
  \label{stress-energy-plane}
\end{equation}
平面应力问题的物理方程为
\begin{align}
  \sigma_x        & =\frac{E}{1-\nu^2}(\epsilon_x+\nu\epsilon_\varphi) \\
  \sigma_\varphi  & =\frac{E}{1-\nu^2}(\epsilon_\varphi+\nu\epsilon_x) \\
  \tau_{x\varphi} & =\frac{E}{1+\nu}\gamma_{x\varphi}
\end{align}
将其代如式\ref{stress-energy-plane}可得：
\begin{equation}
  U=\frac{ER}{2(1-\nu^2)}\int_0^L\int_0^{\rm 2\pi}\int_{-h/2}^{h/2}
  \left(\epsilon_x^2+\epsilon_\varphi^2+2\nu\epsilon_x\epsilon_\varphi
  +\frac{1-\nu}{2}\gamma_{x\varphi}^2\right)
  {\rm d}z{\rm d}\varphi{\rm d}x
  \label{stress-energy-plane2}
\end{equation}

根据Donnell–Mushtari理论\upcite{lee2015free}，应变可以表示为
\begin{align}
  \epsilon_x       & =\frac{\partial \xi}{\partial x}-z\frac{\partial^2\omega}{\partial x^2}                                                                           \\
  \epsilon_\varphi & =\frac{1}{R}\frac{\partial \eta}{\partial \varphi}+\frac{\omega}{R}-\frac{z}{R^2}\frac{\partial^2\omega}{\partial \varphi^2}                      \\
  \tau_{x\varphi}  & =\frac{\partial\eta}{\partial x}+\frac{1}{R}\frac{\partial \xi}{\partial \varphi}-\frac{2z}{R}\frac{\partial^2\omega}{\partial x\partial \varphi}
\end{align}
将上式代入\ref{stress-energy-plane2}即可求出系统的应变能。

\subsection{刚度矩阵及质量矩阵}
根据拉格朗日法，系统自由振动的动力学方程可写为
\begin{equation}
  \frac{\rm d}{{\rm d}t}\frac{\partial T}{\partial \dot{\psi}_i}+\frac{\partial U}{\partial \psi_i}=0
\end{equation}
将$i$取遍全部广义坐标即可得到系统的动力学方程。
考虑到线性系统的动力学方程可以表达为矩阵形式，且$\frac{\rm d}{{\rm d}t}\frac{\partial T}{\partial \dot{\psi}_i}$对应惯性项；
$\frac{\partial U}{\partial \psi_i}$对应刚度项，因此必然有
\begin{equation}
  \frac{\rm d}{{\rm d}t}\frac{\partial T}{\partial \dot{\psi}_i}=\sum_{j=1}^{n}m_{ij}\ddot{\psi_j}
\end{equation}
其中,$m_{ij}$表示质量矩阵第$i$行第$j$列的元素。

同理，对于刚度矩阵也有
\begin{equation}
  \frac{\partial U}{\partial \psi_i}=\sum_{j=1}^{n}k_{ij}\psi_j
\end{equation}
其中$k_{ij}$表示刚度矩阵第$i$行第$j$列的元素。

\subsection{假设振型及边界条件}
由于位移边界条件为固定约束，固柱壳端面的位移及转角均为0。
易知$\frac{\partial \omega}{\partial x}$可表示以$\phi$为轴的旋转；$\frac{\partial \xi}{\partial x}$可表示以$x$为轴的旋转。因此在固定端\songkuohao{$x=0$}，应有
\begin{align}
  \xi                                & =0 \\
  \eta                               & =0 \\
  \omega                             & =0 \\
  \frac{\partial \omega}{\partial x} & =0 \\
  \frac{\partial \xi}{\partial x}    & =0
\end{align}
由于$\xi$分量只要求函数值为0，因此采用一端固支杆的振型作为其$x$方向振型函数，即
\begin{align}
  \sin( & \frac{{\rm \pi} x}{2L})  \\
  \sin( & \frac{3{\rm \pi} x}{2L})
\end{align}
由于$\eta$及$\omega$分量的函数值及一阶导数均要求为0，故选用振型函数为
\begin{align}
  x^2 \\
  x^3 \\
  x^4
\end{align}


在$\varphi$方向上，根据有限单元法的结果，柱壳呈节径型振动。由于各阶振型在周向上呈现周期性，因此应使用三角函数作为假设模态。

文献\cite{lee2015free}中推荐使用$\cos n \varphi$表示$\xi$及$\omega$分量的$2n$节径振动；
使用$\sin n \varphi$表示$\eta$分量的$2n$节径振动。本文中对有限单元法的计算结果进行检查，
发现在各阶振型中，$\xi$及$\omega$分量振动的波峰刚好与$\eta$分量振动的波谷重合，
因此文献中给出的周相位移模式是合理的。
从结构的对称性角度分析，对于子午面，$\xi$及$\omega$是对称分量；而$\eta$是反对称分量，所以这种模态假设方式也是有依据的。

有限单元解中，前三阶振型的节径数分别为6、8、4；
周期分别为$\frac{\rm 2\pi}{3}$、$\frac{\rm 2\pi}{4}$、$\frac{\rm 2\pi}{2}$。
因此，本文采用$\cos(2\varphi)$、$\cos(3\varphi)$、$\cos(4\varphi)$作为$\xi$及$\omega$分量在$\varphi$方向的振型函数；
采用$\sin(2\varphi)$、$\sin(3\varphi)$、$\sin(4\varphi)$作为$\eta$分量在$\varphi$方向的振型函数。

在实际计算中发现，不同节径数的模态是互不耦合的，即在任意一阶模态中，节径数不同的振型函数只有一个非零。
因此，本文对$\omega$分量选用以下振型函数
\begin{align}
  \cos(n\varphi)x^m
\end{align}
其中$n= \rm 2,3,4$,  $m=2,3,4$。
$\eta$分量选用以下振型函数
\begin{equation}
  \sin(n\varphi)x^m
\end{equation}
其中\textit{n}=2,3,4，$m=2,3,4$。
对$\xi$分量选用以下振型函数
\begin{equation}
  \cos(n\varphi)\sin( \frac{m{\rm \pi} x}{2L})
\end{equation}
其中$n=2,3,4$，$m=1,2$。

将$n$依次取3、4、2以计算前三阶振型。

\chapter{计算结果}
前三阶互异振型的频率见表\ref{freq-mod-ana}；其振型见图\ref{level1-ana}-图\ref{level3-ana}。
\begin{align}
  1+1=2
\end{align}

\begin{table}[htbp]
  \centering
  \caption{假设模态解的前三阶互异模态频率}\label{freq-mod-ana}
  \renewcommand\arraystretch{2}
  \begin{tabular}{ccc}
    \toprule
    学科 & \makecell[c]{成绩     \\（百分制）} & 排名 \\\midrule
    \makecell[c]{语           \\文} & 85 & 3 \\
    数学 & 90              & 2 \\
    \bottomrule
  \end{tabular}
\end{table}

\begin{table}[htbp]
  \centering
  \caption{假设模态解的前三阶互异模态频率}\label{freq-mod-ana2}
  \begin{tabular}{lll}
    \toprule
    阶数 & 频率\songkuohao{Hz} & 相对有限元解的误差\songkuohao{\per{}} \\\midrule
    1  & 1567.8 & 5.1           \\
    2  & 1746.6 & 6.1           \\
    3  & 2403.9 & 3.4           \\\bottomrule
  \end{tabular}
\end{table}

\begin{figure}[htbp]
  \centering
  \includegraphics[height=150pt]{./images/mod1-ana.png}
  \caption{假设模态法一阶模态}\label{level1-ana}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[height=150pt]{./images/mod2-ana.png}
  \caption{假设模态法二阶模态}\label{level2-ana}
\end{figure}

\begin{figure}[htbp]
  \centering
  \includegraphics[height=150pt]{./images/mod3-ana.png}
  \caption{假设模态法三阶模态}\label{level3-ana}
\end{figure}
\section{总结}
通过对比假设模态法和有限单元法的计算结果，可以看出，假设模态法对频率的估计与有限单元法比较接近，且误差随振动节径数增加而增大；二者的振型一致。
振动频率的估计值总是偏大，这可能是因为假设模态法强制振型为有限个给定的函数，因此模型受到了额外的“约束”，导致振动频率增加。
}

\backmatter{
  \bibliographystyle{GBT7714--2015}
  \bibliography{ref}
}

\end{document}